Porównać w języku Julia reprezentację bitową liczby 1/3 dla Float16, Float32, Float64 oraz liczby,
która jest inicjalizowana jako Float16, a potem rzutowana na Float64.
a = Float16(1/3)
b = Float32(1/3)
c = Float64(1/3)
d = Float64(Float16(1/3))
println("Float16(1/3) = ", a, " ", bitstring(a))
println("Float32(1/3) = ", b, " ", bitstring(b))
println("Float64(1/3) = ", c, " ", bitstring(c))
println("Float64(Float16(1/3)) = ", d, " ", bitstring(d))
Float16(1/3) = 0.3333 0011010101010101 Float32(1/3) = 0.33333334 00111110101010101010101010101011 Float64(1/3) = 0.3333333333333333 0011111111010101010101010101010101010101010101010101010101010101 Float64(Float16(1/3)) = 0.333251953125 0011111111010101010101000000000000000000000000000000000000000000
Zbadać, jak zmienia się odległość między kolejnymi liczbami zminnoprzecinkowymi reprezentowanymi w komputerze za pomocą języka Julia. Narysować wykres używając Plots zależności odległości od wartości liczby dla zakresu od 1.0 do 1000000.0.
using Plots
x = 1:100:1000000
dist32(x) = eps(Float32(x))
dist64(x) = eps(Float64(x))
y32 = dist32.(x)
y64 = dist64.(x)
scatter(x, [y32 y64], xlabel="x", ylabel="distance", title=["Float32(x)" "Float64(x)"], legend=false, layout=2, colour=[:red :blue])
Jedną z bibliotek numerycznych, jaką dodatkowo będziemy używać na zajęciach jest GSL (język C). Opis jak używać . Korzystając ze wsparcia dla wyświetlania reprezentacji liczb zmiennoprzecinkowych zobaczyć jak zmienia się cecha i mantysa dla coraz mniejszych liczb. Zaobserwować, kiedy matysa przestaje być znormalizowana i dlaczego?
Kod załączyć jako komórka Markdown sformatowana jako C (link). Wynik także jako Markdown (kod albo fragment zrzutu ekranu).
#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>
int main (void)
{
float x = 1e-32;
while(x > 0) {
x = x/(2.0);
printf ("x = ");
gsl_ieee_printf_float(&x);
printf ("\n");
}
return 0;
}
x = 1.10011111011000100011111*2^-108
x = 1.10011111011000100011111*2^-109
x = 1.10011111011000100011111*2^-110
x = 1.10011111011000100011111*2^-111
x = 1.10011111011000100011111*2^-112
x = 1.10011111011000100011111*2^-113
x = 1.10011111011000100011111*2^-114
x = 1.10011111011000100011111*2^-115
x = 1.10011111011000100011111*2^-116
x = 1.10011111011000100011111*2^-117
x = 1.10011111011000100011111*2^-118
x = 1.10011111011000100011111*2^-119
x = 1.10011111011000100011111*2^-120
x = 1.10011111011000100011111*2^-121
x = 1.10011111011000100011111*2^-122
x = 1.10011111011000100011111*2^-123
x = 1.10011111011000100011111*2^-124
x = 1.10011111011000100011111*2^-125
x = 1.10011111011000100011111*2^-126
x = 0.11001111101100010010000*2^-126
x = 0.01100111110110001001000*2^-126
x = 0.00110011111011000100100*2^-126
x = 0.00011001111101100010010*2^-126
x = 0.00001100111110110001001*2^-126
x = 0.00000110011111011000100*2^-126
x = 0.00000011001111101100010*2^-126
x = 0.00000001100111110110001*2^-126
x = 0.00000000110011111011000*2^-126
x = 0.00000000011001111101100*2^-126
x = 0.00000000001100111110110*2^-126
x = 0.00000000000110011111011*2^-126
x = 0.00000000000011001111110*2^-126
x = 0.00000000000001100111111*2^-126
x = 0.00000000000000110100000*2^-126
x = 0.00000000000000011010000*2^-126
x = 0.00000000000000001101000*2^-126
x = 0.00000000000000000110100*2^-126
x = 0.00000000000000000011010*2^-126
x = 0.00000000000000000001101*2^-126
x = 0.00000000000000000000110*2^-126
x = 0.00000000000000000000011*2^-126
x = 0.00000000000000000000010*2^-126
x = 0.00000000000000000000001*2^-126
x = 0
// mantysa przestaje być znormalizowana po przekroczeniu zakresu cechy
Na przykładzie wybranego algorytmu niestabilnego numerycznie:
Nie wolno pokazywać przykładów z wykładu (lub bardzo podobnych)!
Wszystkie punkty przedstawić w postaci notatnika Julii.
# Funkcja:
f(x) = cos(x) ^ 2
# Algorytm niestabilny numerycznie:
# Funkcja obliczająca pochodną funkcji
function f_derivative(f, x, e)
return (f(x) - f(x-e)) / e
end
# 1:
println("Oczekiwany wynik: -2 * cos(1) * sin(1) ~= ")
println(-2 * cos(1) * sin(1))
println()
println("Wyniki algorytmu niestabilnego:")
println(f_derivative(f, 1, 0.001))
println(f_derivative(f, 1, 0.0000001))
println(f_derivative(f, 1, 0.0000000000001))
println(f_derivative(f, 1, 0.00000000000001))
println(f_derivative(f, 1, 0.000000000000001))
println(f_derivative(f, 1, 0.0000000000000001))
println()
# 2. Dobranie niewłaściwej wartości epsilon (e) powoduje zwiększanie błędu.
# Zbyt duże e, rzędu ~0.001 powoduje odbiegający od rzeczywistości wynik algorytmu, a zbyt małe e mniejsze od systemowego epsilon, może spowodować "wyzerowanie" wyniku.
# 3. Wersja stabilna:
function f_derivative2(f, x)
e = sqrt(eps(Float64))
return (f(x+e) - f(x-e)) / (2e)
end
println("Wynik algorytmu stabilnego:")
println(f_derivative2(f, 1))
Oczekiwany wynik: -2 * cos(1) * sin(1) ~= -0.9092974268256818 Wyniki algorytmu niestabilnego: -0.9097129673252824 -0.909297467877046 -0.9092726571680032 -0.8992806499463768 -0.8326672684688673 0.0 Wynik algorytmu stabilnego: -0.9092974290251732